AD-774  723 


NONLINEAR  OPTIMIZATION  USING  THE 
GENERALIZED  REDUCED  GRADIENT  METHOD 

Leon  S.  Lasdon,  et  al 

Case  Western  Reserve  University 


Prepared  for: 

Office  of  Naval  Research 

October  1973 


DISTRIBUTED  BY: 


National  Technical  Information  Service 
U.  S.  DEPARTMENT  OF  COMMERCE 

5285  Port  Royal  Road,  Springfield  Va.  22151 


ikmLm 


o 


DOCUMENT  CONTROL  DATA  -  R  8.  t> 

< Seet.rtry  cJmttificMiion  o(  Ilf/*,  body  ot  mbttr+ot  end  indexing  enno'Atlon  muei  o*  entered  when  the  overmli  report  is  c ietttUed) 


i! 


•Cinatjnc  activity  (Corporate  author) 

Case  Western  Reserve  University 


2*.  REPORT  SECURITY  C  L  ASSl  PfC  A  T*©N 

Unclassified 


2b.  CROUP 


1  REPORT  TITLE 

"Nonlinear  Optimization  Using  the  Generalized  Reduced  Gradient  Method" 

A  DESCRIPTIVE  NOTES  (Type  of  report  etui  incht&lve  dele  a) 

Technical  Report,  1973 

s  authorisi  (Fleet  nesr.e,  middle  tnltimt ,  leet  neme) 

Leon  S.  Lasdon,  Richard  L.  Fox,  Margery  W 

Ratner 

5  REPORT  DATE 

V«.  TOT  AC  MO.  or  PACES 

76  MO.  OF  REFS 

7 

$/t.  contractor  jrant  no. 

9*.  ORIGINATOR'S  REPORT  NUMBER!*)  1 

N00Q14-67-A-0404-0010 

b.  PROJECT  NO 

c. 

96.  OTHER  REPORT  NO(S|  (Any  other  number*  thet  eney  be  eielfred 
thie  report) 

d. 

Technical  Memorandum  No.  325 

10  DISTRIBUTION  STATEMENT 

Unlimited  Distribution 

II  SUPPLEMENTARY  NOTES 

12.  SPONSORING  MILITARY  ACTIVITY 

None 

Office  of  Naval  Research 

13.  A05TRAC  T 


*» 


Generalized  Reduced  Gradient^  methods  are  algorithms  for 
solving  nonlinear  programs  of  general  structure.  This  paper 
discusses  the  basic  principles  of  GRG,  and  constructs  a  specific 
GRG  algorithm.  The  logic  of  a  computer  program  implementing  this 
algorithm  is  presented  by  means  of  flow  charts  and  discussion.  A 
numerical  example  is  given  to  illustrate  the  functioning  of  this 
program. 


Kepro'lucod  by 

NATIONAL  TECHNICAI 
INFORMATION  SERVICE 

U  S  Department  of  Commerce 
Springfield  VA  ??m 


ir ' 


.n 


14  7H 


Nonlinear  Optimization  Using  the  Generalized 
Reduced  Gradient  Method 

by 

Leon  S.  Lasdon 
Richard  L.  Fox 
Margery  W.  Ratner 

Technical  Memorandum  No.  325 


r  ' 
C. 


October  1973 


This  report  was  prepared  as  part  of  the  activities  of  the  Department 
of  Operations  Research,  School  of  Management  and  the  Solid  Mechanics, 
Structures,  and  Mechanical  Design  Division,  School  of  Engineering,  Case 
Western  Reserve  University  (under  Contract  fPWL4-67~A-0404-0010  with  the 
Office  of  Naval  Research).  Reproduction  in  wlfole  or  part  is  permitted  for 
any  purpose  by  the  United  States  Government. 


“asann 


* 

i 


2,  Basic  SSeas  of  GRG 


The  wnlinear  program  to  be  solved  is  assumed  to  have  the  form 
/minimize  f(X)  (1) 

subject  to  (X)  «  0,  i  =  1 . o  (2) 

.  atad  i  Xi  -  ui*  1  “  i; - (3)  ; 

where  X  is  n-vector  and  are  given  lower  and  upper  bounds  u^  >  J 

♦ 

We  assume  m  <  n  since,  in  cost  cases,  m  >_  n  implies  an  infeasible  problem 9 

(1)  -  (3)  * 

or  one  with  a  unique  solution.  .  The  for®  is  completely  general, 


since  inequality  constraints  aav  always  be  transformed  to  equalities,  as  in  (2), 
by  the  addition  of  slack  variables.  The  vector  X  contains  as  components  both 
the  "natural"  variables  of  the  problem  and  these  slacks. 

The  fundamental  idea  of  GRC  is  to  use  the  equalities  (2)  to  express 
m  of  the  variables,  called  basic  variables,  in  terms  of  the  re«aining  n-m  nonbasic 
variables.  This  is  also  the  way  the  Simplex  Method  of  linear  programming 
operates.  Let  X  be  a  feasible  point  and  let  y  be  the  vector  of  basic  variables 
and  x  the  nonbasic  at  X,  so  that  X  is  partitioned  as 

X.  =  (y,x) ,  X  =  (y,x)  (A) 

and  the  equalities  (3)  can  be  written 

g(y,x)  =  0  (5) 

where 

%  =  (gx»  •  •  -  (6) 

Assume  that  the  objective  f  and  constraint  functions  are  differentiable. 
Than, by  the  implicit  function  theorem,  in  order  that  the  equations  (5)  ha  a 
a  solution  y(x)  for  all  x  in  some  neighborhood  of  x,  it  is  sufficient  that 
the  n  x  m  basis  matrix  3g/3y,  evaluated  at  Xt  be  nonsingular . 
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Assume  that  it  is.  Then  the  objective  may  be  expressed  as  a  function 
of  x  only:  < 

F(x)  =  f(y(x),x)  *  (7) 

and  the  nonlinear  program  is  transformed,  at  least  for  x  close  to  x,  to 
a  reduced  problem  with  only  upper  and  Icfwer  bounc.s: 


minimize  F(x) 


(8) 


(9) 


subject  to  » 

^NB  —  X  —  UNB 

/} 

where  and  u^  are  the  vectors  of  bounds  for  x.  GRG  solves  the  original 
problem  (l)-(3)  by  solving  a  sequence  of  problems  of  the  form  (8)-(9). 


Such  problems  may  be  solved  by  simple  modifications  of  unconstrained  minimi- 

t 

zation  algorithms. 

f 

For  the  reduced  problem  (8) -(9)  to  yield  useful  results,  it  is  necessary 
that  x  be  free  to  vary  about  the  current  point  x.  Of  course,  the  bounds 
(9)  restrict  x„  but  it  is  easy  to  move  x  in  directions  which  keep  these 
bounds  satisfied.  The  bounds  on  the  basic  variables,  however,  pose  a  more 
serious  problem.  If  some  components  of  y  are  at  their  bounds,  then  even  a 
slight  change  in  x  L row  x  may  cause  some  bound  to  be  violated.  To  guarantee 
that  this  cannot  happen,  and  to  insure  the  existence  of  the  function  y(x) , 
we  assume  that  the  following  condition  holds: 

Nondageneracy  Assumption 

At  any  point  X  satisfying  (2)-(3),  there  exists  a  partition  of  X  into 
m  basic  variables  y  and  n-ra  nonbasic  variables  x  such  that 

4  <  y  <  «b  (10) 

where  l ^  and  u^  are  the  vector  of  bounds  y  and 

8  ~  ^gpy  is  nonsingular  (IX) 

x’'*’.'!  ass  unset*  °n  is  quite  mild  as  we  show  Icter. 


Consider  now  starting  from  the  feasible  point  X  with  basic  variables 
y  and  nonbasic  variables  and  attempting  to  solve  the  reduced  problem 
(8)-(9).  By  (7)}to  evaluate  the  objective  F(x),  we  must  know  the  values 
of  the  basic  variables  y(x) .  Of  course,  except  for  linear  and  a  few 
nonlinear  cases,  the  function  y(x)  cannot  be  determined  in  closed  form. 
However,  y(x)  can  be  computed  for  any  given  x  by  an  iterative  method  which 
solves  the  equalities  (5) .  Hence  a  procedure  for  solving  the  reduced 
problem  starting  from  X,  is 

(0)  set  i  =  0 

(1)  Substitute  x^  into  (5)  and  determine  the  corresponding  values 
for  y,y^,  by  an  iterative  method  for  solving  nonlinear  eq^c^ions 

(2)  Determine  a  direction  of  motion,  d^  for  the  nonbasic  variables  x 

(3)  Choose  a  step  size  such  that 

*i+i  ■  *1  +  “i  di 

This  is  often  done  by  solving  the  one  dimensional  search  problem 
minimize  F(x^  +  o,di) 

with  a  restricted  such  that  x.^  +  adi  satisfies  the  bounds  on  x,  This 

one  dimensional  search  will  require  repeated  applications  of  step  (1) 

to  evaluate  F  for  various  a  values. 

(4)  Test  the  currant  point  X.^  «  (y^x,)  for  optimal!  :y.  If  not 
optimal,  set  i  =  i  +  \  and  return  to  (1). 

If,  in  step  (1),  the  value  cf  one  or  more  components  cf  y^  exceed 
their  bo\inds,  the  iterative  procedure  must  be.  interrupted.  Fcr  simplicity, 
assume  only  one  basic  variable  violates  a  bound.  Tnen  this  variable  must 
be  mad*.  nonbasic  and  some  component  of  x  which  is  not  on  a  bound  is  made  basi 
After  ♦•h'.s  change  of  basis,  we  have  a  new  function  y(x) ,  a  new  function  F(x), 
and  a  new  redureo  problem.  These  ideas  are  illustrated  geometrically  in 
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Pigure  2.1.  The  initial  point  X  is  on  the  curve  gj  (X)  =0.  We  have 
taken  the  basic  variables  as  (x^x^),  although  the  only  variable  that 
cannot  be  basic  is  x^,  since  it  is  at  lower  bound  of  zero.  The  objective 
of  the  first  reduced  problem  is  F^^^),  which  is  just  the  objective  f  as 
measured  on  the  curve  g2  ~  0  It  possible  that  the  algorithm  minimizing 
F^  might  release  x^  from  its  lower  bound  of  zero,  in  which  case  we  would 
move  interior  to  g2  *  0.  Assume,  for  purposes  of  illustration,  that  this 

dees  not  happen.  Then  we  move  along  g2  *  0  as  indicated  by  the  arrow  until 

we  hit  the  curve  =  0.  At  this  point  the  slack  for  goes  to 

zero.  Since  it  is  basic,  it  must  leave  the  basis,  to  be  replaced  by  one 
of  the  nonbasics,  x2  or  x^.  Since  x^  is  zero,  x,  becomes  basic.  Now  we 
have  a  new  opjective,  F^x^jX^),  with  x^  and  x^  at  lower  bound  of  zero. 

The  algorithm  optimizing  F2  will  determine  that,  if  either  x^  or  x^  is 
released  from  its  lower  bound,  F2  can  be  decreased.  Assume  is  released 
from  its  bound  (actually  and  x^  might  both  be  released  from  their  bounds) . 

Then  the  algorithm  will  begin  to  minimize  F2,  which  is  simply  f  as  measured 

along  the  curve  =  0.  Motion  is  towards  the  x2  axis.  Upon  reaching  it,  x^ 
becomes  zero,  and  another  basis  change  occurs,  with  x-,  becoming  nonbasic 
and  becoming  basic.  Then  optimization  of  the  new  function  F^f-c^x^) 
will  terminate  at  the  constrained  optimum  of  f. 

The  Reduced  Gradient 

- v - 

GRG  can  be  implemented  without  using  derivatives  of  f  or  the  . 

This  requires  methods  for  solving  nonlinear  equations  and  for  minimizing 
nonlinear  functions  subject  to  bounds  which  do  not  use  derivatives.  Although 
such  methods  exist,  there  are  a  much  greater  variety  which  do  require 
derivatives.  The  efficiency  of  these  is  better  understood,  and  their  use 
ir  lar'-:e  nrobl^ns  is  hptter  established. 
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Hence  wa  concern  ourselves  from  now  on  with  GRG  algorithms  which  require 
first  derivatives  of  f  and  g. 

In  minic-i  ing  F  using  derivatives,  we  must  have  a  formula  for  VF t 
F  is  guaranteed  to  be  differentiable  if  f  and  g  are,  and  if  3g/3y  is  non¬ 
singular,  since  then  the  implicit  function  y(x)  is  differentiable.  By(7) 
9F/3xi  =  3f /3xi  +<3f/3y)T  3y/3Xi  (12) 

To  evaluate  3y/3x^,  use  the  fact  that,  if 
gj(y(x),x)  "0,  j  *  l,...,m 
for  all  x  in  some  neighborhood  of  x,  then 

dg^Ah^  =  0  =  (3gj/dy)T3y/3xi  +  3g  /ax^  i  =  l,...,m 
or,  in  matrix  form 

^3g/3y)  3y/3xi  +  3g/3x;L  -  0 
Since  (3g/3y)  is  nonsingular  at  X 

3y/3x±  =  -(3g/3y)_1  3g/3xi  =  B~13g/3xi  (13) 

Using  (13)  in  (12) 

3F/3xi  =  3f/3xi  -  (3f/3y)T  b"1  3g/3xi  (14) 


Let 

x  =  (3f/3y)T  B_1  (15) 

As  is  shown  later,  the  a— vector  t  is  the  Kuhn— Tucker  multiplier 

vector  for  the  constraints  g.  Using  (15),  the  components  of  VF  are 

T 

••  3f/3xi  -  x  3g/3xi  (16) 

Equation  (16)  reduces  to  the  formula  for  the  relative  cost  factors 
in  linear  programming  [1]  if  f  and  all  gj^  are  linear.  Then,  ?f/3x^  =  Cj , 

3f/3y  «  Cg  (the  objective  coefficients  of  the  basic  variables)  and  3g/3x^  «  P^, 
the  column  of  constraint  coefficients  for  x^.  The  vector  u  is  the  simplex 


multiDlier  vector. 


Relation  of  Reduced  Gradient  Formula  and  Kuhn -Tucker  Conditions 

If  X  is  optimal  for  (l)-(3) ,  and  if  the  gradients  of  all  binding 
constraints  of  X  are  independent  (see  [2]),  then  the  Kuhn-Tucker  conditions 
hold  at  X.  To  write  these,  let  u  be  a  Lagrange  multiplier  vector  for  the 
equalities  (2),  and  a  and  3  be  multipliers  for  the  lower  and  upper  bound 
constraints  respectively.  The  Lagrangian  for  (l)-(3)  is 
L  =  f  +  rg  +  a  (/-X)  +  S  (X-ir) 

The  Kuhn-Tucker  conditions^ written  in  terms  of  y  and  x,  are 


3L/3y  *  3f/3y  +  ttb  -  +  8^  ■  0 

(17) 

3L/3x  =  3f/3x  +  tt  3g/3x  “  ax  *  ^x  “  ® 

(18) 

ot  >  0,  8  >  0 

(19) 

a(£~x)  =  8(x-tt)  =  0 

(20) 

where  a^,  3y  are  3ubvectors  of  a  and  6  corresponding  to  the  basic  variables 
y,  and  similarly  for  ax,  If  X  is  optimal,  there  exist  vectors  v, 

a,  6  which,  together  with  X,  satisfy  (17)-(20).  Since  y  is  strictly 
between  its  bounds,  (20)  implies 
a  B  S  “  0 

y  y 

Then  (17)  implies 

tt  =>  -  9f/3y  B  ^ 

so  the  vector  11  in  (15)  is  the  multiplier  vector  for  the  equalities  (2). 

Then  (18)  may  be  written 

3f/3x  +  Tt  3g/3x  -  -  8x  (21) 

The  left  hand  side  of  (21)  is  simply  the  reduced  gradient,  7F(x) . 

To  relate  (21)  to  the  problem  (8) — (9) ,  if  x^  is  strictly  between  its  bounds 

then  a  =3  =  0  by  (20) ,  so 

Xi  xi 

*F/3x,  -  0 


(2?) 
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If  x,  at  lower  bound  then  8  -  0  so 

i  x^ 


?F/3x. 


i  =  o  >0 

Xi" 


(23) 


while  if  x,.  is  at  upper  bound,  -  0  so 

?F/3x.  =  -8  £  0  (24) 

Xi 

But  (22)  -  (24)  are  just  the  optimality  conditions  for  Lhs  reduced  problem 

(8)  -  (9) .  Hence  the  Kuhn-Tucker  conditions  for  (1)  -  (3)  may  be 

viewed  as  optimality  conditions  for  the  reduced  problem  (8)  -  (9) , 

and  ir  in  the  formula  for  the  reduced  gradient  is  the  Kuhn-Tucker  multiplier 

vector.  This  vector  is  useful  for  sensitivity  analysis,  and  GRG  provides 

it  as  a  by-product  of  its  computations. 


Itelatlon  of  Nondegeneracy  Assumption  and  Luenberger 


Constraint  Qualification 

Let  X°  be  an  optimal  solution  to  (1)  -  (3) .  Luenberger 

[2J  has  shown  that  a  sufficient  condition  for  the  Kuhn-Tucker  conditions 
(> 

to  hold  at  X  is  that  the  gradients  of  all  binding  constraints  be  linearly 
independent . 

Assume  that  this  is  the  case.  Then,  at  most  n  of  the  2n  +  ra 
constraints  (1)  -  (3)  can  be  binding  at  X°.  Since  all  m  of  the  equalities 
(2)  are  binding,  at  most  n-m  of  the  constraints  (3)  can  be  binding,  i.e  at 
most  n-m  of  the  variables  X^  can  be  at  a  bound.  Hence  there  will  be  at 
least  m  variables  X.,  satisfying^  <  X°  <  u . .  Consider  now  the  Jacobian 
matrix  of  the  binding  constraints  evaluated  by  X°.  If  all  variables  not  at 
bounds  are  grouped  together,  this  Jacobian  has  the  structure 


m  equality 
tows 

n  -  k 
bound 
rows 


k  >  m  variables 
not  on  bounds 


n  -  k  variables 
on  bounds 


n-k 


J 


Since  this  matrix  must  have  all  ra+(n  -  k)  rows  linearly  independent,  it 

must  have  this  same  number  of  independent  columns.  Since  the  n-k 

rightmost  columns  are  independent,  the  submatrix  J°  must  contain  m 

independent  columns.  Let  B  be  a  nonsingular  m  x  m  submatrix  chosen  from 

J-p  and  let  y  be  the  m-vector  of  ■  .riables  associated  with  the  columns  of  B°, 

with  x  the  vector  cf  the  remaining  n-m  variables.  Then  (  <  v°  <  u 

'  B  J  B 

,  _c  .  , 

vu  3  -s  ncnsingular. 
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That  is,  the  nondegeneracy  assumption  stated  earlier  is  true  at  X°,  so 
it  is  implied  by  Luenbergers  constraint  qualification.  This  information 
is  useful,  since  Luenbergers  qualification  appears  to  be  satisfied  at  the 
optimum,  indeed  at  all  feasible  points,  of  all  but  a  few  pathological 
nonlinear  programs.  However,  problems  can  arise  where  the  binding  constraint 
gradients  become  nearly  dependent,  and  then  B  becomes  nearly  singular, 
and  its  inversion  and  other  operations  with  it  become  numerically  unstable. 

A  computer  program  implementing  GRG  must  test  for  this  near-singularity 
and  attempt  to  correct  it  if  it  occurs. 
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3-  A  GRG  Algorithm 

In  this  section  we  decribe  the  GRG  algorithm  developed  during  the 
period  Nov.  1972-Nov.  1973.  The  major  differences  between  this  algorithm 
and  the  procedure  described  by  Abadie  in  [3]  and  [4 3  are: 

1.  The,  algorithm  works  only  with  the  currently  active  constraints. 
Since,  in  most  problems,  not  all  constraints  are  active,  this 
can  ease  computations  considerably.  The  basis  matrix  has  a  row 
for  each  active,  constraint,  and  changes  size  as  constraints  are 
encountered  or  dropped.  Gradients  of  inactive  constraints  are 
not  required,  a  significant  advantage  in  problems  with  many 
constraints. 

2.  The  algorithm  used  to  optimize  the  objective  on  each  constraint 
intersection  is  the  Davidon-Fletcher-Powell  (DFP'  method  [3],  modi¬ 
fied  to  account  for  upper  and  lower  bounds.  This  should  yield 
more  rapid  convergence  than  the  gradient  or  conjugate  gradient 
procedures  used  by  Abadie. 

3.  A  new  procedure  has  been  constructed  for  deciding  whether 

to  incorporate  a  constraint  into  the  current  constraint  basis. 

The  constraint  is  incorporated  if  the  one-dimensional 
minimum  currently  being  sought  is  on  the  boundary  of  the  current 
constraint  intersection.  Mechanisms  for  determining  this  effi¬ 
ciently  in  the  context  of  GP.G  have  been  developed. 

4.  A  new  basis  change  procedure  is  used.  In  [2],  Abadie  makes  a 
basis  change  if  a  basic  variable  violates  one  of  its  bounds  during 
the  Newton  iteration.  This  can  lead  to  "falsa"  basis  changes  if 
the  Newton  algorithm  is  not  converging,  or  is  converging  but  not 
monotonicallv .  We  wait  until  the  Newton  algorithm  has  converged, 
then  treat  the  violated  bound  as  a  newly  encoxintered  constraint, 
and  apply  the  procedures  in  (3)  above.  This  insures  that  the 
objective  value  after  a  basis  change  is  lower  than  all  previous 
values  (this  is  not  true  in  Abadies  realization) . 

5.  The  one-dimensional  search,  is  the  core  of  the  algorithm  and  Is 
crucial  to  its  efficiency.  We  have  adopted  the  algorithm  described 
in  Ld  to  operate  within  GRG.  This  prccedxtre  is  the  result  of 

.  many  years  of  development,  and  computational  resuits  uc-ing  it  in  un-  * 

c<^mir*4  idrdj&xAti**  hcv*  taw. 

fte  xmWj.  present  $mi  ii&cma  flat*  c&feste  *f  our  G&S  algorithm*  t&U 

context,  discuss  the.  above  ideas  in  more  detail.  The  algorithm  currently 

requires  a  feasible  starting  point.  Work  during  the  next  year  will  include  n 
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designing  a  phase  I  procedure,  which  will  find  a  feasible  point  or 
determine  that  none  exists. 

Let  X  -  (y,  x)  be  a  feasible  po:nt  for  the  constraints  (2)  -  (3). 
Further,  suppose  that  the  first  constraints  are  equalities,  and 
the  remaining  are  inequalities,  with  =  m.  That  is,  the 

constraints  may  be  written: 

g±(X)  =0,  i  “  1, . . . ,m^ 

g^CX)  _>  0,  i  *  m^  +  I . m 

We  define  the  index  set  of  binding  constraints  at  X  as 
IBC  =  {i|g1(X)»  0} 

The  surface  defined  by  a  set  of  active  constraints  will  be  called  the 
constraint  intersection,  S: 

S  =  {x|gjL(X)  *  0.,  ielBC) 

GRG  moves  from  one  such  surface  to  another  as  new  constraints  become  binding 
and  previously  binding  constraints  become  positive.  Ir-  our  realization  of 
GRG,  while  on  a  particular  constraint  intersection,  we  ignore  the  positi  e 
constraints,  except  for  evaluating  them  at  each  point  to  check  that  thev 
are  still  positive.  The  constraint  basis  at  X  contains  a  row  for  each  index 
in  IBC.  Since  the  slacks  of  binding  constraints  are  zero  (i.e  at  lower 
bounds),  there  will  never  be  any  slacks  in  the  basis.  If  there  are  NB 
binding  constraints,  the  NB  basic  variables  are  all  chosen  from  the  n 
"natural"  variables,  X^,...,Xn,  while  the  n  nonbasics  are  the  remaining 
n  -  NB  natural  variables,  plus  the  NB  slacks  of  the  binding  constraints. 

Use  of  Goldfarb  Variable  rfefric  Algorithm 

Since  GRG  requires  at  least  partial  solution  of  a  number  of  reduced 
problems  of  the  form  (8)  -  (9) ,  each  algorithm  for  solving  such  problems 
leads  to  a  variant  of  GRG.  The  choice  of  algorithm  is  critical,  since  once 
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GRG  approaches  an  optimum,  and  no  more  basis  changes  occur,  its  convergence 

rate  is  that  of  the  algorithm  selected.  It  would  be  foolish  to  use  the 

method  of  steepest  descent,  as  its  convergence  rate  is  at  best  geometric, 

with  a  very  small  convergence  ratio  for  problems  whose  Hessian  at  the 

optimum  is  badly  conditioned  [2].  The  conjugate  gradient  method  used 

by  Abadie  is  superlinearly  convergent  [2] ,  but  many  computational 

experiments  have  shown  it  to  be  considerably  slower  in  practice 
(in  terms  of  number  of  Iterations)  than  methods  of  the  variable 

metric  class  [2].  For  this  reason,  we  have  chosen  the  variable  metric 

method  of  Goldfarb  [5],  simplified  for  the  special  case  cf  bounded 

variables,  to  solve  the  reduced  problems. 

The  flow  chart  entitled  "Main  GRG  Program"  illustrates  our 
adaptation  of  Goldfarbs  algorithm.  The  algorithm  is  very  much  as  decribed 
in  [5],  but  simplified  for  the  special  case  of  bounded  variables  (see  [7]). 

The  flow  chart  is  almost  exactly  as  it  would  be  if  the  only  constraints  present 
were  the  bounds  on  the  nonbasic  variables.  All  the  logic  required  to  deal 
with  the  nonlinear  constraints  (2)  is  in  the  one  dimensional  search 
subroutine  on  page  2  of  the  chart.  The  algorithm  chooses  search  directions 
by  the  formula 

di  ,F  <V 

where  H,  is  an  n  x  n  symmetric  positive  semi- definite  matrix.  This 
matrix  projects  any  vector  onto  the  bounds,  i.e,  for  any  vector  v,  Hv  is  zero 
in  the  ith  position  if  is  at  a  bound.  The  initialization  in  block  1, 
cage  1,  and  the  updating  of  block  3,  page  2  (which  forces  row  and  column  r 
of  H  to  zero  when  xr  hits  a  bound)  guarantee  that  H  always  has  this  property. 

The  algorithm  will  minimize  a  positive  definite  quadratic  objective,  subject 
to  upper  and  lower  bounds,  in  a  finite  number  of  iterations.  If  .at  the  optimum, 
n^  n  of  the  variables  are  at  their  bounds,  then,  once  these  variables  reach 
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their  bounds,  the  optimal  values  of  the  remaining  ones  will  be  found  in  at 
nost  n  -  iterations.  Further,  the  nonzero  rows  and  columns  of  H  form  a 
positive  definite  matrix  which  will  converge  (in  the  case  of  quadratic  objective) 
to  the  inverse  Hessian  of  the  function  of  n  -  n3  variables  formed  by  replacing 
the  n,  variables  at  bounds  by  the  values  of  those  bounds. 

Block  2,  page  1  of  the  flow  chart,  is  performed  as  follows: 

e 

The  Kuhn-Tucker  multiplier  for  the  lower  bound  constraint  on  is  =  9F/9x^, 

and  for  the  upper  bound  **  -9F/3x^.  If  x^  is  at  lower  (upper)  bound  and 
Q  u 

A£  (*i^  *s  negative,  then  F  can  be  decreased  by  increasing  (decreasing)  x^, 
i.e.  by  moving  away  from  the  bound.  In  our  program,  from  all  variables  at 
bounds  whose  multipliers  are  negative,  we  choose  the  variable  with  the 

mul tin lier  of  largest  absolute  value.  If  this  value  is  xarger  than  twice 
j  j  djj  ||  (where  j  j  |  j  indicates  Euclidean  norm) ,  we  allow  this  variable 
to  leave  its  bound  by  setting  the  corresponding  diagonal  element  of  H 
(currently  zero)  to  one.  This  causes  ] j  d^  | |  to  increase.  We  then  test 
the  remaining  bounded  variables  for  negative  multipliers,  and  repeat  the 
above  procedure.  The  test  against  | j  j |  insures  that  we  do  not  leave  a 
constraint  subspace  until  | |  d^  | j  becomes  "small  enough".  This  helps  to 
prevent  zigzagging,  where  we  constantly  leave  and  then  return  to  the  same 
subspace. 

Goldfarbs  algorithm  provides  search  directions  for  the  one 
dimensional  search  subroutine,  in  which  the  variables  of  the  pioblem  are 
assigned  new  values-  This  subroutine  finds  a  first  locaJ  minimum  for  the 
problem 

minimize  F(x  +  ctd) 

where  F  is  the  reduced  objective,  x  the  initial  values  of  the  nonbasic 
variables,  an'1  d  the  search  direction.  The  ’deration  subscript, 

has  beer,  dropped  for  convenience.  The  di-ecclon  d  U  always  a  direction 
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of  descent,  i.e 

dT  VF  (?)  <  0 

The  procedure  starts  with  a  search  for  three  a  values,  A,B,  and 
C,  which  satisfy 

0  <_  A  <  B  <  C 
F  (x  +  Ad)  >  F  (x  +  Bd) 
and 

F  (x  +  Cd)  >  F  (x  +  Bd) 

Then  the  interval  [A ,C]  contains  a  local  minimum  of  F(x  +  ad).  In 
block  18  of  page  1  of  the  one  dimensional  search  flow  chart,  a  quadratic 
is  passed  through  A,  B,  and  C,  with  its  minimum  at  D.  On  page  3  of  the 
flow  chart,  the  points  A,  B,C,  D  are.  used  to  initiate  an  iterative  cubic 
interpolation  process  which  yields  the  final  a  value. 

The  logic  on  the  first  two  pages  of  the  flow  chart  locates  the 
points  A,  Bj  C.  In  doing  this,  the  choice  of  initial  step  si^e,  a^, 

(block  1,  page  1),  is  important.  With  Goldferbs  algorithm  or  other 
variable  metric  methods,  aQ  is  set  equal  to  the  optimal  a  value  from  the 
previous  search  except  when  this  causes  too  large  a  change  in  the  variables. 
The  theoretical  basis  for  this  is  that,  as  a  variable  metric  converges, 
the  optimal  a  values  should  converge  to  1,  the  optimal  step  for  Newton's 
Method.  Hence  the  previous  optimal  step  is  a  good  approximation  to  the 
current  one.  This  must  be  modified  when  the  method  is  restarted,  for  example 
when  a  new  constraint  is  encountered '  or  the  basis  is  changed,  since  then 

an  optimal  step  much  lees  than  unity  is  generally  taken.  Hence,  we  require 

-3 

that  the  change  in  any  nonbasic  variable  larger  than  10  in  absolute  value 
not  exceed  .05  times  its  value,  while  the  change  in  any  variable  smaller 
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than  10  in  absolute,  value  cannot  exceed  0.1.  If  the  largest  a  value 
meeting  these  conditions  is  a\  then,  at  iteration  i,  is  given  by 
=  min 

Ihe  loop  2  -  3  -  4  -  5  halves  the  step  size  until  a  value  FB  <  FA 

is  achieved,  or  until  SLPFLG  =  0.  The  variable  SLPFLG  is  initialized 

at  zero  in  subroutine  NEWG  and  is  set  to  unity  in  block  3  of  the  NEWS 

flow  chart  if  a  new  constraint  has  been  encountered  and  the  minimum  along 

d  is  interior  to  that  constraint.  The  test  in  block  6  of  the  one 

dimensional  search  flow  chart  is  false  only  if  the  step  size  has  been 

halved  at  least  once  in2-3-4-5,  in  which  case  K1  is  the  function 

value  corresponding  to  C.  The  test  in  block  7  prevents  the  subroutine 

from  trying  to  pass  a  quadratic  through  3  points  which  are  too  widely 

separated  in  function  value.  It  also  insures  that  the  subroutine  will 

cut  back  the  step  size  if  a  large  function  valve  is  returned  by  block 

3.  This  is  used  to  force  a  cutback  when  the  NEWTON  algorithm  in  block  3 

30 

does  not  converge,  by  setting  FB  to  10  .  The  test  on  FC  in  block  12  has 

the  same  purpose.  If  K1  is  too  large,  then  block  8  gererates  a  new  C 

point  1/3  of  the  distance  from  B  to  C.  Then  the  loop  8-9-10-11-12 

is  traversed  until  FC  is  not  too  Jarge. 

With  R  =  0  (which  occurs  if  and  only  if  (a)  a  K1  or  FC  which  was 

too  large  has  never  been  generated,  or  (b)  SLPFLG  **  0,  the  loop  10  -  14 

transforms  the  points,  A,  B,  C  in  figure  3..1(a)  into  those  shown  in  figure 

3.1(b).  The  step  size  is  doubled  each  time  until  the  points  A,  B,  C 

bracket  the  minimum.  If  K1  or  FC  ever  becomes  too  large,  or  if  SLPFLG  m  1, 

R  is  set  to  1  (block  9).  Then  (10)  -  (14)  transforms  points  as  shown 

In  figures  3.2  (a)  thru  3.2(c).  Instead  of  doubling  the  step,  a  constant 

30 

B  -  A,  is  added.  Since  FC  may  have  been  set  to  10 


rr 
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wher.  NEWTON  did  not  converge  in  block  11,  there  must  be  a  provision 
for  resetting  R  to  0  when  3  steps  of  size  (C-B)/3  have  been  taken. 

This  is  accomplished  by  blocks  15,  16,  and  17. 

The  quadratic  interpolation  in  block  18  yields  a  fourth  point, 

B,  with  function  value  FM,  somewhere  between  A  ar.d  C.  In  block  19,  a 
cubic  polynomial  is  passed  through  the  4  points  FA,  FB,  FC,  FM,  and 
its  minimum  is  located  at  the  point  E.  The  optimality  tests  in  block 

20  are  passed  if  the  percent  difference  becween  (a)  the  F  values  at  the 

current  and  previoue  interpolated  points  and  (b)  the  values  of  F  and 

-4 

the  cubic  at  E  are  sufficiently  small.  Currently  z  -  10  .  In  block 

21  the  usual  situation  is  as  shown  in  figures  3.3  (a)  and  3.3  (b) . 

Removal  of  an  end  point  leaves  4  points  which  bracket  the  minimum  and 
these  are  used  in  the  next  cubic  *it»  If  a  bracket  cannot  be  formed  by 
removal  of  an  end  point,  the  highest  end  point  is  discarded  and  a  new 
cubic  fit  is  performed.  Such  a  situation  is  shown  in  figure  3.4. 

If  the  optimality  test  in  block  20  Is  passed,  and  the  point  E 
is  less  than  a,  then  E  is  accepted  as  optimal.  Otherwise,  F  is  evaluated 
at  a  in  block  22.  If  its  value  there  is  smaller  than  F(x),  a  is  returned. 

If  not,  and  FC  >  FB,  the  quadratic  interpolation  block  is  entered;  otherwise, 
the  subroutine  terminates  with  an  error  message. 

The  blocks  in  the  one  dimensional  search  flow  chart  labeled  "E  valuate 
F  (x  +  ad)"  are  where  the  basic  variables  are  changed  in  response  to 
changes  in  the  nonbasics.  As  shown  in  Che  accompanying  flow  chart, 
this  block  contains  2  subroutines.  Subroutine  REDOBJ  uses  a  Newton  iteration 
to  find  new  values  for  the  basic  variables  y,  given  the  values  x  +  ad  for 
the  nonbasics.  It  also  checks  for  violations  of  currently  nonbinding 
constrain-.''  an^  of  bounds  on  basic  variables,  using  subroutine  3SCHNG 


to  deal  w't”  these  bounds. 


Figure  3.5  -  Tamgent  Vector  Estimate 
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Subroutine  NEWG  decides  whether  or  not  to  stay  on  a  newly  encountered 
constraint.  We  now  proceed  to  discuss  these  subroutines. 

Subroutine  REDOB J  (REDuced  OBJective)  is  the  first  routine  discussed 

thus  far  that  is  not  very  similar  to  a  procedure  for  unconstrained  minimization. 
Its  input  is  a  vector  of  nonbasic  variables  x  =  x  +  ad,  where  a  is  the 

current  trial  value  of  the  step  size,  and  x  is  the  vector  of  nonbasic  variable', 
at  the  start  of  the  one  dimensional  search.  It  solves  the  system  of  NB 
nonlinear  equations 

(y  x  +  ad)  *  0,  i 6  IBC 

for  the  NB  basic  variables  y.  As  in  [3]  and  [4],  this  is  done  using  the 
pseudo-Newton  algorithm 

yt+J  =  yt  “  B_1  <yt>  *  +  «d),  t  -  0,1 . 

where_gg  is  the  vector  of  binding  constraints. 

-1 

The  algorithm  is  called  pseudo-Newton  because  B  is  not  re-evaluated  at  each 
step  of  the  algorithm,  as  in  the  standard  Newton  method.  When  r  is  given 

A 

its  first  trial  value  in  a  particular  one  dimensional  search,  X  is  equal  to 
1 

the  feasible  point  X  with  which  we  began  the  search.  As  long  as  Newton 
A  _ 

converges,  X  remains  equal  to  X,  at  least  until  the  search  is  over.  If 

A 

Newtcn  fails  to  converge,  X  may  be  set  equal  to  the  most  recent  feasible 
point,  and  B  ^  is  recomputed. 

To  explain  blocks  1  and  2  on  page  1  of  the  REDOBJ  flow  chart,  consider 
the  tangent  plane  to  the  constraint  surface  at  X.  This  is  the  set  of  all 
vectors  (a,b)  satisfying 


(3g/8y)a  +  (3g/3x)  b  »  0 
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where  all  partial  derivative  matrices  are  evaluated  at  X.  In  GRG,  the 
change  in  x,  b,  is  given  by 
b  =  ad 

The  corresponding  vector  a  is  called  the  tangent  vector,  v.  Since  any 
scale  factor  multiplying  v  is  unimportant,  we  may  as  well  take  a  *  1, 
yielding 

v  =  (3g/3y)  1  (3g/3x)d  (25) 

In  our  program,  v  is  computed  at  X,  the  initial  point  of  the  one 
dimensional  search.  This  vector  is  used  to  find  initial  values,  y  ,  by 
the  formula 

y  =  y  +  a,  v 

as  illustrated  in  figure  3.5.  Using  these  initial  values,  Newton  finds 
the  feasible  point  X.^,  Then,  at  X^,  v  is  not  recomputed.  The  old  v  is 
used,  but  emanating  now  from  X^,  to  yield  the  next  set  of  initial  values 
as 

yo  =  yx  +  (a2  -  ax)v 

Using  these,  Newton  finds  the  point  Xj  of  figure  3.5.  This  procedure 

is  repeated  until  Newton  fails  to  converge  (or  until  the  one  dimensional 

sear  :h  is  over),  whereupon  v  is  recomputed  at  the  last  feasible  point. 

-1 

Both  this  logic  and  the  logic  of  computing  B  have  the  objective 
of  computing  derivatives  oe  f  and  the  binding  g  ,  and  of  inverting  B,  only 
when  absolutely  necessary.  If  Newton  converges  at  each  point  of  a  one 
dimensional  search,  then  no  derivatives  or  matrix  inversions  are  required 
during  the  search. 

Newton  is  considered  to  have  converged  it  the  condition 

NORMG  -  max  |gt(X  )  !  <  EPSNKWT 
iclBC 

is  met  ITLTN  iterations.  Currently  (using  single  precision 
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arithmetic;  8  decimal  digits),  EPSNEWT  *  10  and  ITLIM  =10.  If 
NORKG  has  not  decreased  in  any  set  of  5  consecutive  iterations  (or  the 
above  condition  is  not  met  in  10  iterations)  Newton  has  not  converged, 
and  the  2  alternatives  on  page  1  are  tried,  in  an  effort  to  achieve 
convergence . 

The  first  alternative  is  tried  if  the  gradients  of  the  objective  and 
the  binding  constraints  have  not  yet  been  computed  at  the  last  feasible 
point,  XPREV.  These  gradients  are  evaluated,  and  an  approximation  to  the 
true  B  ^at  XPREV  is  computed,  using  the  current  B  ^  (evaluated  at  some 
earlier  feasible  point)  and  the  partial  derivatives  evaluated  at  XPREV. 
This  approximate  inverse  is  computed  as  follows: 

-1 

Let  B  be  the  basis  matrix  evaluated  at  XPREV,  and  B  the  basis  inverse 

’  o 

evaluated  at  some  other  point.  The  approximation  to  3  ^  is  based  on  the 
identity 


B  =  B  +  (B  -  B  ) 
o  o 

=  Bq  (I  -  B~*  (Bq  -  B))  (26) 

Let 

A  =  (Bq  -  B)  (27) 

Then,  taking  the  inverse  of  both  sides  of  (25)  yields 

B-1  -  (I  -  A)"1  B~J  (28) 

If  the  points  at  which  Bq  and  B  are  evaluated  are  sufficiently  close 
together,  the  norm  of  A  will  be  less  than  unity,  and  (I-A)  ^  can  be 
expanded  in  a  power  series: 

(I-A)-1  »  I  +  A  +  A2  + .  (29) 

This  serieo  ran  b  i  used  to  approximate,  the  Newton  correction 


6  =  B'1G 


(30) 


vVre  G  is  the  vector  cf  binding  constraints.  Using  (28)  and  (29)  in 


=  [I  +  A  +  A 


.ufV, 


' yields 


th 

The  i  order  approximation  to  5 , 8 ^  is  obtained  by  truncating  the  series 

expansion  above  at  the  term  A*: 

$  ±  =  T.  A?  B^G  ,  i  =  0,1,2,... 

j=o  ° 

The  vectors  may  be  determined  recursively  as  follows: 

$  =  B_1G 

o  o 

=  (i  +  a)So  =SQ  +  a  8  0 

$  2  “  (I  +  A  +  A2)J  o  -  £o  +  A(I  +A)^o  +  a51 

In  general 

5j+l‘*o  +  A*j»  j  0,3., ...  ,i 

or  using  the  definition  of  A  in  (27) 

$i+1  -  SQ  +^j  ~  K1  sSr  j  =  (31) 

Returning  to  alternative  1  on  page  2  of  the  flow  chart,  this  alternative 
is  implemented  by  choosing  the  order  of  approximation  i  (  i  >_  1)  and,  within 
the  Newton  subroutine,  approximating  5  <l8  cL  using  the  recursion  (31). 

For  i  *  1,  this  is  the  approximation  suggested  by  Abadie  in  [3]  -  [4). 

If,  after  trying  alternative  1,  Newton  again  fails  to  converge,  alternative 

2  is  tried.  This  alternative  computes  the  true  B  *  at  XPREV,  uses  it  in 

(25)  to  compute  a  new  tangent  vector,  and  returns  to  the  Newton  subroutine.  If 

Newton  still  fails  to  converge,  the  final  alternative  is  tried:  the  objective 

30 

is  set  to  a  very  large  value  (10  ),  and  we  leave  subroutine  REDOBJ,  returning 

to  the  one  dimensional  search  subroutine.  The  large  objective  value  will  cause 
the  search  subroutine  to  decrease  c.  This  will  occur  either  in  block  3 
or  block  11  of  the  subroutine.  Tills  cutback  procedure  will  continue  until 
Newton  converges.  This  must  nappen  eventually,  since  Newton  will  converge 
if  the  initial  point  is  close  enough  to  the  solution. 


-22- 


Once  Newton  has  converged,  we  check  the  positive  constraints  to 
sea  if  any  are  now  binding  or  violated.  Let  the  point  Newton  has  obtained 
be 


where 


A  A  _  a 

x  =  (y(a),  x  +  cd) 


gj^yCcb.  X  +  ad)  •-  0,  ielBC 


Assume  that  some  nonbinding  constraints  are  violated; 
gi(X)  <  0,  ielVC 

The  program  attempts  to  find  the  point  at  which  the  first  nonbinding 
constraint  went  to  zero.  That  is,  we  wish  to  find  the  smallest  value  of 
a,  a*,  such  that  all  constraints  in  IBC  are  binding,  exactly  one 
constraint  from  IVC  is  binding,  and  all  other  constraints  in  IVC  are 
positive.  Hence  a*  satisfies 

gA  (y(«),  x  +  ad)  -0,  ieXBC  (32) 

gk  (y(a),  x  +  ad)  «  0 
and 

(y(ct),  x  +  ad)  >  0,  ielVC,  i  i  k 

where  keIVC.  This  is  illustrated  in  figure  3.6  where  IVC  *  (2,3),  k  »  2. 
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Cf  course,  the  index  k  is  not  known  in  advance,  s..  linear  interpolation 

is  used  in  block  3,  page  2,  to  estimate  the  values  o£  a  *  and  k.  The 
Jacobian  for  (32),  J,  is 

J  ~  8  )  S 

l  ^  tf 

where 

w  »  9gk/3y 
t  =  (3gk/3x)Td 

T 

and  s  is  an  N3  component  column  vector  whose  elements  are  (3g^/3x)  d 
for  i  £ IBC.  Since  J  involves  only  adding  a  border  to  the  current  basis 
matrix,  B,  its  inverse  is  easily  computed  if  B  ^  is  known.  In  our  program, 
the  border  vectors  w,s,t  are  evaluated  at  the  last  feasible  point,  XPREV, 
and,  as  a  first  try,  the  current  B  ^  is  used,  even  if  it  was  evaluated  at 
a  point  other  than  XPREV.  The  resulting  J  1  may  be  a  kind  of  “hybrid”, 
but  if  Newton  converges,  an  inversion  of  B  has  been  saved.  Looking  at 

A  —1 

figure  3.6,  since  Newton  converged  to  X,  using  the  current  B  and  starting 
from  XPREV,  one  would  expect  it  to  converge  to  the  point  X*,  which  is  closer 
to  XPREV,  ever,  though  an  additional  constraint  has  been  added.  If  Newton 
fails  to  converge,  the  same  three  alternatives  as  before  (with  minor 
modifications  -  see  block  5,  page  3)  are  tried,  until  convergence  is  achieved. 
Then  the  remaining  constraints  in  IVC  are  checked  to  see  if  they  are  positive. 
If  one  or  more  are  negative,  then  the  linear  interpolation  estimate  of  which 
constraint  was  violated  first  was  in  error.  We  go  back  to  the  linear  inter¬ 
polation  block,  as  the  first  step  toward  finding  a  new  value  of  a*.  This 
cycle  may  be  repeated  a  number  of  times,  but  the  sequence  of  a*  values 
should  decrease.  If  not,  the  procedure  Is  not  working  and  we  stop  with  an 
error  message.  Otherwise,  we  leave  this  section  of  code  with  k  as  the  index 
cf  r  new  binding  constraint  and  «*  ns  the  a  value  at  which  this  constraint 
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is  binding. 

The  last  major  task  in  RED03J  is  to  deal  with  any  basic  variables 
which  have  violated  their  bounds.  This  occurs  in  subroutine  BSCHNG 
(block  4,  page  2  of  REDOB J  flow  chart).  Turning  to  its  flow  chart,  the 
first  action  (block  1*  page  1)  is  to  check  if  any  bounds  on  basic 
variables  have  been  violated.  Note  that,  if  any  g^  constraints  had  been 
violated,  the  basic  variables  are  now  equal  to  their  values  at  the 
point  X*  =  (y*,  x  +  a*d) .  For  example,  in  figure  3.6,  the  y  value  at 
X*  might  be  below  its  lower  bound.  If  some  basic  variables  do  violate 
their  bounds,  we  proceed  through  essentially  the  same  logic  as  is  used 
in  dealing  with  violated  constraints.  The  result  is  a  value  of  a, 
such  that  all  components  of  y^)  except  one  are  strictly  between 
their  bounds,  with  that  one,  Xjj,  equal  to  one  of  its  bound  values. 

Then,  after  storing  the  current  feasible  point  as  XPREV,  we  leave 
subroutine  REDOBJ. 

The  next  subroutine  encountered  in  evaluating  the  reduced  objective 
is  NEWG.  If  a  new  g^  constraint  and/or  a  bound  on  z  basic  variable  has 
been  made  binding  in  REDOBJ,  NEWG  decides  whether  or  not  it  should  remain 
binding.  The  behavior  of  the  objective  F(x)  is  the  determining  factor. 
Using  figure  3.6  as  an  example,  if  the  one  dimensional  minimum  of  F 
occurs  for  a  less  than  a*,  then  the  new  constraint  g^  is  not  made  binding. 
Tho  one  dimensional  minimum  lies  in  the  current  constraint  intersection, 
so  the  step  size  is  reduced,  and  the  search  for  an  optimal  a  continues. 

If,  however,  F  is  still  decreasing  at  a*,  the  one  dimensional  minimum 
is  on  the  boundary  of  the  current  constraint  intersection.  The  new 
constraint  is  made  binding  and  the  one  dimensional  search  terminates. 
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If  both  a  bound  on  a  basic  variable  and  a  new  constraint  become  binding 
in  REDOB J,  this  logic  still  applies,  but  with  a*  replaced  by  min  (a*,  c^). 

The  program  has  been  designed  to  make  these  decisions  without  using 
derivatives.  If  in  block  1,  page  1,  the  objective  at  min  (a*,  c^)  is 
larger  than  that  at  the  last  feasible  point,  then  the  nev.  constraint  is  not 
aduw.  This  is  case  1  of  figure  3.7.  If  the  objective  is  smaller,  then 
we  must  determine  if  case  2a  or  case  2b  of  figure  3.7  holds.  The  reduced 
objective  is  evaluated  at  a  point  nine-tenths  of  the  distance  between  the 
last  feasible  o  value  and  min  (a*,  ce^).  If  the  objective  value  there  is 
smaller  than  that  at  min  (a*,  a^),  case  2b  is  assumed  to  hold  and  the  new 
constraint  is  not  added.  Otherwise  case  2  a  holds.  Either  a  new 
constraint  is  added  or  the  basis  is  changed,  after  which  we  return  to  the 

start  of  the  main  GRG  program  with  a  new  reduced  problem  to  solve. 

The  new  constraint  incorporation  or  basis  change  is  carried  out  in 

subroutine  CONSBS.  The  input  to  this  subroutine  is  a  list  of  indices  of 
variables,  the  candidate  list.  In  block  2  of  the  NEWG  flow  ch'.rt,  this 
list  is  set  equal  to  the  current  list  of  basic  variables.  Then  CONSBS 
is  called.  The  outputs  of  CONSBS  are  (a)  a  new  list  of  binding  constraint 
indices  (b)  a  new  list  of  basic  variable  indices  and  (c)  a  new  basis  inverse > 
called  BIKV  in  the  flow  chart  of  CONSBS.  On  page  1  of  this  flow  chart, 
the  array  IREM  contains  the  list  of  rows  which  remain  to  be  pivoted  in. 

This  is  initialized  in  block  1.  The  subroutine  operates  in  2  modes,  indicated 
by  the  variable  MODE.  When  MODE  »  1,  CONSBS  will  choose  pivot  columns 
from  whatever  candidate  list  was  input  to  it.  If  a  basis  inverse  could 
not  be  constructed  from  columns  in  this  candidate  list,  or  if  the  original 
candidate  list  included  all  variables  (NCAND  «*  N,  block  2),  MODE  is  set 
to  2.  and  '"ONSB^  *;ill  choose  pivot  columns  from  the  list  of  all  admissable 
columns.  A  column  is  admissable  if  its  variable  is  farther  than  EPSBOUNDS 
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( currently  10  )  from  its  nearest  bound,  and  if  it  has  not  yet  been  pivoted 

in. 

The  main  loop  of  COSSBS  begins  at  block  3.  A  pivot  row  is  chosen  as 
IROW  in  block  4.  The  choice  of  ISV  in  block  5  is  motivated  by  the  desire 
to  have  basic  variables  as  far  from  their  bounds  as  possible,  so  that  fewer 
basis  changes  will  be  required.  The  other  criterion  influencing  the  choice 
of  basic  vari'Aixes  is  that  the  basis  matrix  should  be  well-conditioned. 

We  try  to  insure  this  by  choosing  as  a  prospective  pivot  column  that  index, 
I,  in  ISV  yielding 

max  |TAB  (IROW,  I) | 
icISV 

Tills  is  done  in  block  6.  If  the  element  chosen  passes  2  tests  we  pivot 

on  it  (block  8) ,  transforming  the  Jacobian  and  entering  that  column  into 

B  ■*".  The  column  pivoted  in  is  marked  inadmissable  (block  7),  and  the 

■~1 

procedure  is  repeated  for  each  binding  constraint  until  either  B  has  been 
constructed  (N  branch,  block  9)  or  the  candidate  list  has  been  exhausted 
(Y  branch,  block  10). 

The  two  tests  that  a  pivot  element  must  pass  are  (a)  its  absolute 
value  must  be  larger  than  EPSPIVOT  (currently  10_^)  and  (b)  the  absolute 
value  of  the  ratio  of  all  other  elements  in  the  pivot  column  to  the 
pivot  element  must  be  less  than  RTOL  (currently  100).  The  first  test 
insures  that  we  do  not  pivot  on  an  element  that  is  essentially  zero,  while 
the  second  protects  against  the  generation  of  elements  of  large  absolute 
value  in  B  1.  Such  values  are  symptomatic  of  ill-conditioned  basis  matrices. 
If  either  test  is  failed  while  MODE  =  1,  we  simply  do  not  pivot  ir.  the 
current  row,  and  move  on  to  the  next  one. 
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£f,  when  node  1  terminates,  3  has  not  yet  been  constructed,  we 
attempt  to  complete  its  construction  by  considering  columns  not  in  the 
original  candidate  list.  Mode  2  is  entered  at  marker|5|of  the  flow  chart. 

V 

In  block  11,  the  candidate  list,  ICAND,  is  reset  to  the  set  of  all 
remaining  adtnissable  columns,  MODE  is  set  to  2,  and  we  return  to  the  start 
of  the  main  iterative  loop.  If,  in  this  second  phase,  a  pivot  element  fails 
the  absolute  value  test,  we  temporarily  mark  all  columns  in  ISV 
inadmissable  (by  setting  their  indicator  in  the  IGNORE  array 
to  1,  in  block  12)  and  choose  a  new  ISV  array.  If  all  adnissable  matrix 
elements  in  a  row  fail  the  absolute  value  test,  the  matrix  is  considered 
to  be  singular.  If  a  pivot  element  fails  the  ratio  test  in  mode  2,  it 
is  deleted  from  ISV  in  block  13,  and  the  ISV  element  with  second  largest 
absolute  value  is  tried,  and  so  on.  If  all  fail,  we  set  ISKIP  *  1  in  block 
14,  which  causes  us  to  pivot  on  the  element  originally  selected,  ignoring 
the  ratio  test. 
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4.  A  Numerical  Example 

Consider  the  problem 

minimize  F(X)  =  (X^  -  1)  +  (X^  -  0.8) 

subiect  to 

gl(X)  =  xx  -  X2  -  X3  o 
g2(X)  =  -x2  +  X2  -X4  *  o 

g3(X)  =  Xx  +  X2  -X5  -  1  *  0 

X^o,  o  <  X2  <  0.8,  x3  0,  X4  >  0,  X5  >_  0 
where  X^,  X4,  X^  are  slack  variables.  The  feasible  region  for  this 
problem  is  graphed  in  figure  4.1,  along  with  contours  of  constant  value 
of  the  objective.  The  constrained  optimum  is  on  the  surface  G2  -  0  at  the 
point  X^  =  0.894,  X2  =  0.8.  The  starting  point  is  X^  =  0.6,  X2  =  0.4, 
which  lies  on  the  line  =  0,  and  X2  is  the  initial  basic  variable. 

Hence 

y  =  X2,  x  =  (x1,x2)  =  (Xp  X5),  B  =1 
The  objective  of  the  first  reduced  problem,  F  (x) ,  is  obtained  by  solving 
g3  for  X2,  yielding 

x2  -  1  +  x5  -  xx 

and  substituting  this  into  f,  yielding 

F  (x)  =  (Xx  -  l)2  +  (0.2  +  X5  -  Xx)2 
whose  gradient  at  x  =  (0.6,  0)  is 

7F(x)  =  (3F/3x1  ,  3F/3x2)  =  (0,  -  0.8) 

Since  is  at  lower  bound,  the  initial  H  matrix  is  (see  block  1,  page 
1,  main  GRG  program  flow  chart) 


d  -  -  H  “  F  =  (0,0) 
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In  block  2,  page  1  of  the  main  GRG  program,  ve  check  if  X^  should  be 

released  from  its  lower  bound.  Since  the  Kuhn-Tucker  multiplier  for 

this  bound  is  negative: 

/ 

=  3F/3x2  =  -  0.8 

x2  =  Xj  is  released  from  its  lower  bound.  The  new  H  is 


so  now 

d  =*H  ?  F  -  (0,  0.8) 

and  we  enter  the  one-dimensional  search  subroutine  with  cT  <*  +®. 

This  initial  one  dimensional  search  varies  X1  and  X5  according  to 

Xx  «  0.6 

At-  *=  0  +  0.8a 

so  we  move  along  the  vertical  line  shown  in  figure  4,1.  The  initial 

step  size  is  chosen  to  limit  the  percent  change  in  any  variable  to  less 

than  5%  or,  if  a  variable  is  zero,  to  limit  the  absolute  change  in 

such  variables  to  be  less  than  0.1.  This  latter  criterion  applies 

here,  since  only  X2  is  changing,  and  its  initial  value  is  zero.  Hence, 

the  initial  a,  a  ,  is  chosen  such  that 
o’ 

0.8a  »  0.1 

o 

or 

a  -  0.125 
o 

The  sequence  of  a  and  objective  values  generated  by  the  one  dimensional 
search  is: 


—1 

a 

— 

cbiective 

0 

0.32 

0.125 

0.25 

0.250 

0.20 
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X 


The  seauence  of  points  generated  is  shown  as  points  1,  2,  3  in  figure  { 

\ 

4.1,  For  a  =  0.250,  =  0.6,  which  lies  on  the  constraint  =0. 

The  fact  that  g1  has  become  binding  is  detected  in  block  6,  page  2  of 
RED03J.  Since  the  reduced  objective  is  lower  at  a  =  0,250  than  at 
the  previous  point,  the  N  branch  is  taken  in  block  1  of  NEWG.  The 
reduced  objective  is  computed  at  a  *  0.2375,  and  since  its  value  there, 

0.2041, is  larger  than  the  value  at  a  *  0.250,  the  one  dimensional 
minimum  over  the  interval  [0,  0.250]  is  assumed  to  occur  at  a  =  0.250, 
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This  corresponds  to  motion  along  the  line  g^  *  0.  The  sequence  of 
a  and  objective  values  generated  is 


objective 

0 

0.2 

0.025 

0.1658 

0.05 

0.1352 

0.10 

0.0848 

0.20 

X2  <  0.8 

violated 

Figure  4.1  shows  the  corresponding  values  of  and  X2  as  points 
5  thr*  8.  The  last  a  value,  a  *  0.20,  yields  X^,  =  0.84,  X2  =  0.84, 
which  violates  the  upper  bound  on  X2«  This  is  detected  in  REDOBJ  by 
subroutine  BSCHNG,  which  attempts  to  satisfy  the  upper  bound  by 
solving  the  system 

g^  (y,  x  +  ad)  =  0.6  +  1.2a  ~  X2  =  0 

X2  =  0.8 

which  yields  a  =  0.166,  corresponding  to  point  9  of  figure  4.1.  The 
objective  value  at  point  9  is  0.0400,  which  is  smaller  than  the  value  of 
0.848  at  the  last  feasible  point,  point  7.  Hence,  in  NEWG,  we  take 
the  N  branch  in  block  I,  and  evaluate  F  at  a  =  0.160  in  the  next  block, 

A 

(point  10,  figure  4.1)  yielding  F  =  F  =  0.0433.  This  is  larger  than 
the  value  at  point  9,  so  9  is  accepted  as  the  minimum,  and  subroutine 
C0NSBS  is  called.  This  yields  a  new  set  of  basic  variables,  and  a  new 
basis  as  follows: 

IBC  -  {1} 

basic  variables  y  = 
nonbasic  variables  x  =  (X2,  X^) 

B_1  =  (1) 

After  leaving  C0NSBS,  we  return  to  page  1  of  the  main  GRG  flow  chart 
to  begin  the  third  major  iteration . 

^ince  both  nonbasic  variables  are  at  upper  bounds,  H  is  set 

7,oro  matrix.  To  obtain  the  current  reduced  objective,  we  solve 
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the  binding  constraint  for  tie  basic  variable  in  terms  of  the  nonbasic: 

*1  (30  -  Xt  -  %  -  X3  =  0 

SO 

-  X2  +  X3 

Substituting  the  above  into  the  objective  yields  the  reduced  objective  as 
F(x)  =  (J^  +  X3  -  l)2  +  (X,  -  0.8) 2 
whose  gradient  at  X2  =  0.8,  **  0  is 

VF  <=  <-0.4,  -  0.4) 

In  block  2,  page  1  of  the  main  GRG  program,  X^  has  a  negative  multiplier 
with  value  -0.4,  so  it  is  released  from  its  lower  bound. 


and  the  search  direction  is 


d  =  -H7F  -  (0,  0.4) 

We  begin  the  one  dimensional  search,  with  the  nonbasics  varied 
according  to 

X2  =  0.8 

X3  <*  0  +  (0.4)o 

The  a  and  objective  values  generated  by  this  search  are 


a 

obiective 

0 

0.16 

0.32 

0.04 

0.018 

g2  violated 

The  corresponding  values  of  X^  and  X2  are  shown  as  points  11  and  12 
in  figure  4.1.  Subroutine  RED0B J  detects  the  violation  of  g2»  and 
attempts  to  make  g?  binding  by  solving  the  system 
Cj  (y,  x  +  ctd)  =  X^  -0.8-0. 4a  =0 
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2 

g2  Cy»  x  +  ad)  =  -X^  +  0.8  =  0 

This  is  accomplished  in  one  iteration  of  Newton's  method,  starting 
from  initial  values 

(Xx  ,  a)  =  (0.893,  0.234) 

(computed  in  block  3,  page  2  of  REDOBJ),  with  inverse  Jacobian 


1 -  - 

1  -0.4 

0  -0.578 

-1 

J 

*c 

! 

-1.728  0 

-2.50  -1.44 

The  final  values  are  *  0.894,  a  *  0.236  which  corresponds  to  point 
13  of  figure  4.1.  The  objective  value  at  that  point  is  0.0111,  which 
is  lower  than  the  value  at  the  last  feasible  point,  point  11.  Hence, 
in  NEWG,  the  objective,  is  evaluated  at  point  14,  which  is  0.9  of 
the  distance  between  points  11  and  13,  corresponding  to  a  *  0.2285. 

A 

The  value  there  is  F  »  0.0117,  which  is  larger  than  the  value  at  point 
13,  so  subroutine  CONSBS  is  called  to  construct  a  new  basis.  The 
results  are 

IBC  =  {2} 

basic  variables  y  * 

nonbasic  variables  x  =  (X2,  X^) 

B-1  =  -0.5590 

Returning  to  page  1  of  the  main  GRG  program  flow  chart,  the 
reduced  gradient  is 

VF  -  (-0.118,  0.118) 

Since  X2  is  at  upper  bound  and  at  lower  bound,  the  Kuhn-Tucker 
multipliers  for  these  bounds  are 

Xj  =  =  0.118 

/ 

V2  -  ^F/3X4  =  0.118 


Since  both  are  positive,  point  13  is  optimal. 
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